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Approximate  analytical  solutions  have  been  obtained  for  surface  temperature  and  heat 
flux  for  the  case  of  two  embedded  temperature  sensors.  The  solutions  have  been  verified 
and  the  range  of  validity  has  been  established  using  several  methods  including  comparisons 
with  an  exact  analytical  solution  for  a  linear  problem  and  numerical  calculations  for  a  non¬ 
linear  problem.  The  performance  of  the  model  is  presented  in  both  the  frequency  domain 
and  the  time  domain.  A  propagation  of  error  analysis  is  presented  and  is  used  to  establish 
the  optimum  spacing  between  the  sensors.  The  solutions  place  no  restrictions  on  the 
boundary  or  initial  conditions  and  rely  only  on  current  values  of  temperature  and  its  rate  of 
change.  The  temperature  dependence  of  transport  properties  is  accounted  for  in  an 
approximate  way.  The  method  requires  minimal  computation  and  is  suitable  for 
implementation  as  a  real-time  sensor. 
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Nomenclature 

Fourier  number,  at / x 2 
Thermal  conductivity 
Temperature 

Rate  of  change  of  temperature 

Time 
Fleat  Flux 

Distance  from  heat  flux  surface 
Thermal  diffusivity 
Phase  Angle 
Angular  frequency 


I.  Introduction 

Tide  measurement  of  transient  heat  flux  and  surface  temperature  in  heat-sink  combustion  chambers  continues  to 
present  technical  challenges  to  the  instrumentation  engineer.  Sensor  failure  rates  are  high  and  measurement 
accuracies  and  uncertainties  are  not  well  characterized.  These  shortcomings  have  had  a  significant  impact  on  some 
recent  programs  which  have  used  heat  sink  test  articles  to  acquire  data  for  the  validation  of  heat  transfer  predictions 
at  liquid  rocket  engine  operating  conditions1'2. 

There  are  numerous  types  of  heat  flux  sensors  but  a  relatively  small  subset  is  capable  of  operating  in  rocket 
chamber  conditions  where  heat  flux  levels  can  exceed  108  W/m2  and  surface  temperatures  of  1000  K  are  typical. 
Diller3  reviewed  the  devices  that  have  been  used  and  organized  them  into  methods  that  rely  on  temperature 
differences  over  a  spatial  distance  with  known  thermal  resistance  and  temperature  differences  over  time  with  known 
thermal  capacitance.  The  most  commonly  used  method  has  been  the  coaxial  thermocouple  which  is  an  example  of 
the  second  type.  A  thermocouple  junction  is  formed  on  the  surface  of  the  chamber  between  a  wire  of  one  type  of 
thermocouple  material  and  a  surrounding  sheath  of  another  type.  The  heat  flux  is  determined  from  the  measured 
temperature  boundary  condition  using  a  one-dimensional  transient  solution  to  the  heat  equation.  The  junction  is 
typically  very  thin  and  is  often  formed  by  lightly  scratching  the  surface  to  drag  filaments  of  one  type  of  material 


1  Propulsion  Research  Engineer,  Aerophysics  Branch,  Space  and  Missiles  Propulsion  Division,  Propulsion 
Directorate,  10  East  Saturn  Blvd. 
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across  the  electrically  insulating  layer  to  the  other  type.  In  some  applications,  when  erosion  of  the  surface  occurs  the 
junction  is  continuously  reformed  and  this  has  led  to  the  description  of  coaxial  thermocouples  as  “eroding 
thermocouples.”  However,  in  heat-sink  chambers  it  is  quite  common  to  find  that  the  junction  disappears  at  some 
point  during  a  test  and  the  sensor  fails. 

Other  methods  have  been  developed  which  do  not  rely  on  surface  temperature  measurements  but  embed  the 
sensors  within  the  wall  where  they  are  protected  from  erosion.  In  the  null-point  calorimeter  a  hole  is  drilled  from 
the  backside  of  the  chamber  wall  and  a  thermocouple  is  inserted.  The  bead  is  brazed  or  resistance  welded  to  the 
bottom  of  the  hole.  “Null-point”  refers  to  a  distance  from  the  bottom  of  the  hole  to  the  inner  wall  where  the 
disturbance  to  the  flow  of  heat  caused  by  the  hole  results  in  the  junction  reading  nearly  equal  to  the  inner  wall 
temperature.  The  construction  of  null-point  calorimeters  is  challenging.  The  junction  cannot  be  visually  inspected 
and  large  measurement  errors  can  result  from  manufacturing  flaws4. 

Another  method  using  embedded  temperature  sensors  is  the  plug-type  heat  flux  gauge  of  Liebert5.  An  annular 
groove  is  machined  into  the  chamber  wall  to  form  a  post  and  thermocouples  are  attached  at  several  axial  points 
along  the  outside  of  the  post.  A  polynomial  curve  is  used  to  extrapolate  the  temperatures  to  the  wall  position  and  an 
integral  method  is  used  to  calculate  the  total  heat  load  to  the  plug  from  transient  temperature  measurements.  Two- 
dimensional  effects  can  be  significant  in  this  type  of  device.  The  dimensions  of  the  groove  are  critical  and 
significant  errors  can  result  from  the  disturbance  to  the  flow  of  heat6. 


II.  Derivation  of  Model 


In  the  following  analysis  we  adopt  the  approach  of  Liebert5  of  extrapolating  a  polynomial  to  the  surface.  The 
temperature  profile  in  the  body  (see  figure  1.)  is  approximated  using  a  power  series  in  x  with  time  dependent 
coefficients. 

T(x,t)  =  ^jCi(t)x"  (1) 

1=0 

However,  in  addition  to  matching  the  temperatures  at  the  measurement  points,  the  polynomial  is  also  required  to 
match  the  second  derivatives  with  respect  to  x.  These  values  are  obtained  from  the  time  derivatives  of  temperature, 
which  are  assumed  to  be  available  as  experimental  measurements,  and  the  heat  equation. 


dT  d2T 

- =  a - T 

dt  dx 


(2) 


Adding  these  additional  constraints  allows  a  higher  order  polynomial  to  be  used  and  markedly  improves  the  time 
response  of  the  sensor.  Explicit  expressions  for  the  c:  (t)  can  be  obtained  by  substituting  (1)  into  the  right  side  of 

(2)  then  solving  for  the  cAt)  in  terms  of  the  measured  temperatures  and  their  rates  of  change.  For  two 

measurement  points  there  will  be  four  measured  quantities  and  the  power  series  can  be  carried  to  the  fourth  term.  If 
second  derivatives  are  available,  the  series  can  be  carried  to  the  sixth  term  and  so  on. 


Figure  1  Heat  flux  sensor  geometry 
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Evaluating  the  series  at  the  surface,  x=0,  and  combining  terms  in  order  to  make  the  linear  dependence  on  the 
measurements  apparent,  the  temperature  and  heat  flux  expressions  can  be  written  as  follows, 


T  =T 
1o  1\ 


i. 


-+r. 
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(4) 


In  deriving  these  solutions,  no  assumptions  have  been  made  with  respect  to  the  boundary  conditions  at  the  top  of 
the  block  and  in  fact  it  is  completely  arbitrary.  For  example,  the  top  surface  could  be  actively  cooled.  Also  the 
second  sensor  at  x?  does  not  need  to  be  embedded  within  the  block  and  can  be  located  at  the  top  surface.  The  sensor 
at  x i  may  be  located  at  the  x=0  surface.  Furthermore,  we  have  made  no  assumptions  with  respect  to  the  origin  of 
time.  The  solution  relies  only  on  the  current  values  of  temperature  and  rates  of  change  of  temperature.  Referring 
back  to  the  categorizations  of  methods  defined  by  Diller3,  this  method  is  a  hybrid  of  the  two  types  as  it  relies  on  both 
spatial  and  temporal  variations  in  temperature. 

The  method  used  to  evaluate  the  rate  of  change  of  temperature  is  critical  to  the  success  of  the  technique.  A 
simple  finite  difference  calculation  will  almost  certainly  result  in  an  unacceptable  level  of  noise  in  the  derivative.  A 
polynomial  smoothing  filter,  also  known  as  a  Savitzky-Golay7  3  filter,  can  be  effective  in  reducing  noise  and  can  be 
calculated  very  efficiently.  The  filter  is  implemented  as  a  convolution  in  the  time  domain. 
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T(t.)= 

i=-m 

(5) 

m 

(6) 

i=—m 


Substituting  (5)  and  (6)  into  (3)  and  (4)  we  arrive  at  the  final  forms  for  the  surface  temperature  and  heat  flux 
expressions. 


m 

T0{tn)=  X  eiTl{tn+i)+ 

i  =  -m 

m 

Q  0(O  =  -*(r0)E  SiTx{tn+i)  +  hiT2  (tn+i) 

i=—  m 

The  coefficients  are  defined  as  follows. 
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6«j  (x2  —  Xj) 


(7) 
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h ,  =  a , 


■  + 


x. 


b  2x[  +  2x{x2  -  x 
6 a2(x2  -  x, ) 


(12) 


Note  that  once  the  coefficients  have  been  calculated,  only  4m+2  multiplications  and  2m+l  additions  are  required 
to  evaluate  T0  and  q0 ,  or  about  100  processor  operations  are  required  for  a  typical  value  for  m  of  10. 

The  uncertainties  ST0  and  Sq0  can  be  evaluated  by  treating  each  temperature  measurement  as  an  independent 

random  variable  and  using  the  “square  root  of  the  sum  of  squares”  approach.  We  also  assume  that  uncertainties  in 
thermocouple  locations  and  material  properties  are  negligible. 


,  1/2 


\i=—m 

/  \  1/2 

'  m  \ 

Sq0(t„)  =  k(Tj  Xg/^+Zt,.2^2 

\i=-m  J 


(13) 


(14) 


Minimum  values  of  ST0  and  Sc/(l  can  be  achieved  by  making  an  optimum  choice  for  the  order  of  the  polynomial 
used  to  obtain  the  ai  and  hi .  For  a  fixed  level  of  uncertainty  in  the  measurements,  the  coefficients  of  a  high  order 
method  will  result  in  a  larger  value  for  the  total  uncertainty;  however,  in  the  presence  of  high  frequency  content,  a 
high  order  polynomial  will  fit  the  trend  more  closely  and  reduce  the  values  of  STn ,  thereby  reducing  the  total 

uncertainty.  The  optimum  order  can  be  found  by  increasing  the  order  until  a  minimum  in  the  uncertainty  is  found. 
However  this  increases  the  amount  of  computation  and  in  practice  we  find  that  the  quadratic  smoothing  filter  is  a 
good  compromise  and  can  be  used  as  the  default  method.  The  coefficients  for  the  quadratic  smoothing  filter  are  the 
following. 


fl,  = 


Z-/ 

j=—m  j=-m 


1/ 

j=—m 


2m  + 1 


2j2  £/ 
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(15) 


(16) 


ITT.  Verification  of  Model 

To  illustrate  the  behavior  of  the  approximate  solutions,  we  make  comparisons  with  the  exact  solution  for  the 
temperature  within  a  semi-infinite  slab  exposed  to  an  oscillating  heat  flux  at  the  surface. 
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When  (18)  and  its  time  derivative  are  substituted  into  (3)  and  (4)  the  results  are  linear  combinations  of  sine  and 
cosine  terms  of  the  form: 


T0  =  Acos(cot  +  (j){)  +  Bcos(cot  +  </)2)  +  Csm(cot  +  ^,)  +  Dsm(cot  +  (j)-,)  (19) 


The  four  terms  can  be  combined  into  a  single  term  of  the  following  form. 


T,,  =  Gain 


n  f  \ 


a 

\co) 


Ml 


sin(<y?  +  Phase ) 


(20) 


The  amplitude  and  phase  behaviors  are  shown  in  figures  2  and  3.  For  greater  generality,  the  results  are  plotted  as 
functions  of  two  non-dimensional  numbers.  The  x-axis  is  a  non-dimensional  angular  frequency  which  is  the 
reciprocal  of  a  Fourier  number  based  on  the  angular  frequency  of  the  heat  flux  and  the  distance  from  the  surface  to 
the  first  sensor,  Fo~ 1  =  x\col  a ,  and  the  individual  curves  are  for  various  values  of  the  ratio  x2/xi.  These  two 
variables  completely  determine  the  solutions.  Considering  the  surface  temperature  first,  we  see  immediately  that  the 
solution  has  the  characteristics  of  a  low  pass  filter  but  with  a  region  of  positive  gain  in  the  range  of  angular 
frequencies  from  1  to  20.  We  also  see  that  the  magnitude  of  the  gain  increases  with  the  ratio  x2/x!.  The  phase  of  the 
solution  begins  to  lag  at  angular  frequencies  above  1  except  for  x2/X|~4  where  there  is  a  range  from  0.4-4  where  the 
phase  angle  is  positive. 

The  behavior  of  the  solution  for  heat  flux  is  given  in  figure  3  and  is  substantially  similar.  In  this  case,  the  peak 
in  the  gain  increases  with  x2/x!  but  is  significantly  lower  than  the  peak  in  the  surface  temperature  gain.  In  fact,  for 
heat  flux,  the  gain  behavior  is  quite  a  weak  function  of  x2/x i  over  the  range  from  1.5-4.  The  phase  behavior  is  also  a 
weak  function  of  x2/xi  with  all  values  showing  a  roll  off  beginning  near  an  angular  frequency  of  0.5. 

In  addition  to  correctly  reproducing  the  amplitude  and  phase  of  a  heat  flux  signal,  the  model  must  also  be  robust 
in  the  presence  of  noise.  We  now  consider  the  effect  of  the  parameter  x2/xi  on  the  propagation  of  error  from  noise  in 
the  temperature  measurements.  Equations  (13)  and  (14)  are  expressions  for  the  uncertainties  in  surface  temperature 
and  heat  flux  when  the  smoothing  filter  approach  has  been  used.  An  alternative  approach  for  deriving  these 
expressions  is  to  start  with  (3)  and  (4)  and  assume  that  the  uncertainties  in  temperature  and  rate  of  change  of 
temperature  are  known  quantities.  This  approach  yields  succinct  expressions  that  reveal  the  role  of  x2/X[.  We 
assume  the  error  in  each  measurement  is  independent  and  the  thermocouple  positions  and  material  properties  are 
error  free.  The  sensitivity  coefficients  are  simply  the  coefficients  in  (3)  and  (4).  If  we  assume  5TX  ~  ST2 ,  <5T,  *  S f, 
and  a,  «  a2 ,  the  uncertainty  of  the  surface  temperature  is, 


2  2 


FT0  = 


(x2  -xxy 


-ST2  +■ 


( 5x2  -  8X[X2  +  5*!  ) 

3 6a2 (x2  -Xj)2 


£T: 


The  uncertainty  for  the  heat  flux  is, 


(21) 


$r2  + 


(5x,4  +  4xj3x2  +  4xjX2  +  5x2 )  ^p2 
36cr(x2  -Xj)2 


(22) 
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j  - x2/x1=1.5  - x2/x1  =2.0  -  -  -  x2/x1=3.0  —  -  x2/x1=4.0  | 


- x2/x1=1.5  - x2/x1  =2  -  -  -  x2/x1=3  —  -  x2/x1=4 


Figure  2  Surface  temperature  gain  and  phase 


Equations  (21)  and  (22)  show  that  uncertainty  is  a  function  of  the  positions  of  the  sensors,  and  therefore  we  can 
anticipate  that  there  will  be  optimum  locations  that  minimize  uncertainty.  To  make  quantitative  predictions  we  need 
estimates  for  ST  and  ST  and  these  were  obtained  from  experimental  data  that  will  be  presented  later.  The  specific 
values  used  here  were  5T  =1  K  and  ST  =  20  K/s.  Figure  4  shows  the  uncertainty  in  surface  temperature  based  on 
(21)  and  in  fact  it  has  a  minimum  at  x2/x!=3.6,  but  the  minimum  is  quite  broad  and  any  value  from  2.3  to  5  would 
work  as  well.  This  region  corresponds  well  to  the  region  of  maximum  accuracy  identified  previously.  Flowever,  the 
uncertainty  increases  rapidly  for  x2/xi<2  so  this  range  should  be  avoided.  The  minimum  in  the  heat  flux  uncertainty 
(Figure  5)  based  on  (22)  is  somewhat  sharper  with  the  optimum  value  2.3.  The  heat  flux  uncertainty  has  been 
normalized  by  the  magnitude  of  the  heat  flux  so  it  can  be  seen  that  the  minimum  in  the  uncertainty  is  approximately 
0.75%  of  the  measured  value.  Again,  values  of  x2/xi<2  should  be  avoided.  Based  on  these  uncertainty  estimates 
plus  the  gain  and  phase  results  described  above,  the  optimum  value  for  x2/xj  is  2. 3-2. 5. 
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IV.  Temperature  Dependent  Properties 

If  thermal  properties  are  temperature  dependent  the  heat  equation  is  non-linear  and  analytical  solutions  such  as 
(18)  are  not  possible.  For  this  case  the  model  was  verified  and  validated  using  numerical  results.  We  chose  as  a  test 
case  an  idealized  version  of  the  typical  test  conditions  we  see  in  our  laboratory:  a  step  change  in  heat  flux  from  0- 
4x1 07  W/m2.  The  MatLab  non-linear  PDF  solver,  pdepe  was  used  to  generate  temperature  histories  for  the  surface 
and  internal  points  and  these  were  converted  back  to  heat  flux  and  surface  temperature  using  the  model.  Two  cases 
of  x2/x i  are  shown.  For  both  cases  the  sensor  nearest  the  surface  is  located  at  a  depth  of  1.905  mm  (0.075  inch).  The 
total  thickness  of  the  slab  was  41.91  mm  (1.65  inch)  and  the  material  was  high-conductivity  copper.  The  heat  flow 
was  assumed  to  be  one-dimensional.  The  grid  contained  131  evenly  spaced  points  and  the  temperatures  were  saved 
at  one  millisecond  intervals.  Results  using  a  261  point  grid  differed  by  approximately  0.1  K  and  did  not  change  the 
model  results. 

The  ability  of  the  model  to  reproduce  the  surface  temperature  and  heat  flux  boundary  conditions  is  shown  in 
figure  6.  For  both  cases  of  x2/xi,  after  0.1  seconds  the  polynomial  model  has  converged  to  within  0.1  K  for  the 
surface  temperature  and  within  1%  of  the  heat  flux,  but  the  magnitude  of  the  overshoot  is  significantly  smaller  for 
x2/x!= 2.3  than  3  as  we  should  expect  from  the  gain  plots  presented  earlier. 
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V.  An  Application  of  the  Model 

To  illustrate  the  application  of  the  model  we  present  data  taken  in  a  sub-scale,  heat-sink,  rocket  combustion 
chamber.  Type  K  thermocouples  were  embedded  in  a  high  conductivity  copper  block  that  formed  one  side  of  a  2.54 
cm  square  channel.  The  locations  of  the  sensors  were  those  given  above  in  the  description  of  the  temperature- 
dependent  properties  calculation.  The  reactants  were  gaseous  hydrogen  and  gaseous  oxygen  at  a  mass  flow  rate 
ratio  of  1:5.5.  The  propellants  were  completely  mixed  and  reacted  at  the  measurement  location.  The  pressure  of  the 
chamber  was  42  bar  (600  psig)  and  the  mean  velocity  of  the  gas  was  approximately  250  m/s  (820  ft/s).  In  figure  7 
the  upper  figure  shows  the  temperatures  measured  within  the  block  labeled  as  “Front  Side”  for  xt  and  “Backside” 
for  x2.  Also  shown  are  results  for  surface  temperature  based  on  (3)  and  uncertainty  based  on  (13)  which  is  read  on 
the  right-side  vertical  axis.  The  reference  in  the  legend  to  Tsurfi  and  T3  refers  to  the  data  from  the  third  of  six 
measurement  positions.  In  this  example,  which  used  an  early  version  of  the  technique,  the  uncertainties  of  the 
temperatures  and  rates  of  change  of  temperature  were  obtained  from  the  uncertainties  in  linear  regression  parameters 
over  a  data  window  60  milliseconds  wide  and  containing  61  measurements.  This  is  equivalent  to  a  polynomial 
smoothing  filter  based  on  a  linear  function  in  distinction  to  the  quadratic  filter  recommended  above  and  given  in 
(15)  and  (16).  The  uncertainties  increase  in  the  regions  where  the  rates  of  change  of  temperature  are  not  constant  as 
occurs  at  the  beginning  and  end  of  the  test  and  are  largely  due  to  the  fact  that  the  linear  regression  was  used.  The 
uncertainties  propagate  to  the  surface  temperature  calculation  resulting  in  the  spikes  at  1.6  and  2.2  seconds.  The 
heat  flux  results  are  shown  in  the  lower  figure.  The  maximum  value  reached  was  2.46xl07  W/m2  (15.0  Btu/in2/s). 
The  uncertainty  in  heat  flux  due  to  noise  during  the  middle  of  the  run  is  approximately  0.5%  of  the  measured  value. 
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- x2/x1=3  - x2/x1=2.3 


Figure  6a  Effect  of  temperature  dependent  properties  on  surface  temperature  prediction. 
Test  case  is  step  change  in  heat  flux  from  0  to  4e7  W/m2  (24.5  Btu/in2/s). 


- x2/x1=3  - x2/x1=2.3 


Figure  6b  Effect  of  temperature  dependent  properties  on  accuracy  of  heat  flux.  Test  case 
is  a  step  change  in  heat  flux  from  0-4e7  W/m2  (24.5  Btu/in2/s).  Heat  flux  converges  to 
within  1%  after  0.1  seconds. 


VI.  Conclusion 

A  method  has  been  described  for  measuring  surface  temperature  and  heat  flux  based  on  two  temperature  sensors 
embedded  in  the  wall  of  a  chamber.  The  method  does  not  require  surface  junction  thermocouples  which  are  prone 
to  failure  and  produce  noisy  signals  in  rocket  engine  flows  and  is  well  suited  for  studies  of  the  effects  of  surface 
features  on  heat  transfer  enhancement.  A  polynomial  extrapolation  approach  is  used  where  the  function  is 
constrained  to  match  the  temperatures  and  the  second  derivatives  with  respect  to  distance  at  the  measurement  points. 
The  method  requires  only  current  values  of  temperature  and  its  rate  of  change  and  the  boundary  and  initial 
conditions  are  arbitrary.  The  algorithm  can  be  represented  as  a  low  pass  fdter  and  the  gain  and  phase  behavior  have 
been  characterized.  The  placement  of  the  sensors  affects  the  frequency  cutoff  and  the  noise  response  and  optimum 
values  for  the  relative  positions  of  the  sensors  have  been  obtained.  The  method  utilizes  time  derivatives  of 
temperature  data  which  can  be  efficiently  calculated  using  a  polynomial  smoothing  filter.  The  method  is 
computationally  efficient,  requiring  approximately  100  multiply  and  add  operations  and  is  suitable  for 
implementation  in  a  digital  signal  processor. 
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—  Tsurf3  K  — FrontSide  — BackSide  — T3-Noise 


—  HF3  —  HF3-Noise 


Time  (s) 


Figure  7  Surface  temperature  and  heat  flux  from  experimental  data  obtained  in  a  subscale 
rocket  combustion  chamber 
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